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We consider a normal metal - superconductor (N-S) junction in the regime, when electrons in 
the normal metal are driven out of equilibrium. We show that the non-equilibrium fluctuations 
of the electron density in the N-layer cause the fluctuations of the phase of the order parameter 
in the S-layer. As a result, the density of states in the superconductor deviates from the BCS 
form, most notably the density of states in the gap becomes finite. This effect can be viewed as a 
result of the time reversal symmetry breaking due to the non-equilibrium, and can be described in 
terms of a low energy collective mode of the junction, which couples normal currents in N-layer and 
supercurrents. This mode is analogous to the Schmid-Schon mode. To interpret their measurements 
of the tunneling current, Pothier et. al [Phys. Rev. Lett. 79, 3490 (1997)] had to assume that 
the energy relaxation rate in the normal metal is surprisingly high. The broadening of the BCS 
singularity of the density of states in the S-layer manifest itself similarly to the broadening of the 
distribution function. Mechanism suggested here can be a possible explanation of this experimental 
puzzle. We also propose an independent experiment to test our explanation. 
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I. INTRODUCTION 

In metals the inelastic scattering rate l/ri„ at low 
enough energies is determined by electron-electron inter- 
actions. In a clean Fermi liquid er e > h this inelastic rate 
can be estimated as h/ri n ~ e 2 ftp, where e is the energy 
of the quasiparticle, r e is the elastic scattering time, and 
tp is the Fermi energy. This familiar result reflects only 
the phase volume of the final state for an inelastic pro- 
cess, while the corresponding matrix element is an en- 
ergy independent constant. In the dirty limit er e < h 
the inelastic rate is significantly enhanced as compared 
with the clean case due to long range diffusive correla- 
tions of-jSingle electron wave-functions in the disordered 
systemaa, see Refs. || for more detailed discussion. 

The inelastic scattering rate is not by itself an observ- 
able quantity. However, inelastic collisions of electrons 
profoundly affect the behavior of the system and, there- 
fore, l/rj n in many cases can be extracted from experi- 
mental data. For example, from magnetoresistance one 
can evaluate quantitatively the dephasing time t v , which 
often coincides with l/rj„. The dephasing time describes 
the loss of phase coherence, as electrons move diffusively 
in the bulk of a metallic sample. This loss of coherence 
cuts off otherwise divergent weak-localization correction. 
The dephasing time has been extensively studied, and 
we believe that the existing theory, allows a good under- 
standing of the experimental dataB. 

Another effect of inelastic scattering is the energy re- 
laxation described by the time r e . This is the time it 
takes for a " hot" quasiparticle with energy e much larger 
than temperature T to thermalize with all the other elec- 
trons. Theoretically, it is given bytl 
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(1.1) 



where G m (L) oc L is the dimensionless (in units of 
e 2 /2tt%) conductance of a sample of size L, D is the dif- 
fusion constant and d the dimensionality of the system. 
This result follows from the Fermi Golden Rule, but now 
the phase volume is multiplied by the matrix element, 
which is no longer a constant. Not only this matrix el- 
ement substantially depends on the energy transferred, 
but it even diverges at small energies due to the wave 
function correlation, see Refs. [| 

To determine r e experimentally, one has to apply an 
external perturbation to drive the system out of equi- 
librium, and then to measure the distribution func- 
tion of electrons. Recently, an elegant and important 
experiment^ was performed to measure directly the elec- 
tronic distribution function /(e). An external voltage 
U applied to a copper wire caused an electric current 
J, thus driving the wire out of equilibrium. In order 
to determine the non-equilibrium distribution function 
the authors of Ref. [| fabricated an additional electrode 
connected with the wire by a tunneling contact. The 
distribution function /(e) was extracted from the mea- 
surements of the tunneling conductance Gt(V) of this 
contact as a function of the bias voltage V using the pro- 
cedure as follows. Assuming that the density of electronic 
states in the wire is energy independent, one can present 
the tunneling conductance Gt(V) as the convolution of 
/(e) with the tunneling density of states in the additional 
electrode p(e) 



(1.2) 
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As follows from Eq. (1.2), the more pronounced is the en- 
ergy dependence of the density of states in the additional 
electrode p(e), the more precisely one can determine the 
distribution function /(e) measuring Gt{V). For this 
reason the authors of Ref. ^ used a superconducting elec- 
trode to take advantage of the BCS singularity in the 
density of states p(e) . The existence of this singularity in 
the equilibrium was convincingly determined by indepen- 
dent measurements at J = 0. The data on the tunneling 
conductance Gt(V) were fitted by Eq. (1.2) yielding the 
distribution function /(e) and, thus, the energy relax- 
ation time. 

This procedure produced quite unexpected results. 
First of all, the extracted relaxation time turned o ut t o 
be two orders of magnitude shorter than that of Eq. (1.1). 
Moreover, no dependence of the relaxation time on the 
energy e was observed. This would mean the failure of 



the theory lying behind the derivation of Eq. (1.1). 

However, the theory is based on regular expansion in 
inverse powers of the dimensionless conductance G m (see 
Ref. H and references therein). The observation ota good 
metallic conductance G m S> 1 in the experiment).] justi- 
fies the applicability of this theory. 

In this paper we,attempt to explain the puzzling results 
of the experimentu by lifting the main assumption in the 
interpretation of the data - independence of the density 
of states in the superconductor of the electronic distri- 
bution in the normal metal. We calculate the tunneling 
conductance between the superconducting and metallic 
films explicitly and find that interaction of the tunneling 
electrons with non-equilibrium fluctuations of the cur- 
rent in the normal layer smears the BCS singularity in 
the density of states in the superconductor. This effect 
has nothing to do with energy relaxation. Nevertheless, 
it effectively broadens the energy dep ende nce of the ex- 
perimental /(e), extracted from Eq. fll.2] ) with the use 
of the equilibrium BCS density of states p(e), whereas 
the real distribution function remains sharp. As a result, 
the energy relaxation time appears much shorter, than it 
really is. We found that this effect is not small as inverse 
dimensionless conductance of the normal metal 1 / G m . 
Instead, the magnitude of the effect is proportional to the 
inverse conductance of the superconductor in the normal 
state 1/G S . Under the condition G m ^> G s , which we 
assume in the present paper, this effect dominates the 
real energy relaxation. 

The remainder of the paper is organized as follows. 
In Sec. [n], we present a phenomenological derivation of 
our main results. In the same section we suggest an in- 
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dependent experiment to test our theory. Section 
is devoted to the rigorous analysis of the the tunneling 
density of states under non-equilibrium conditions. Our 
findings are summarized in Conclusions. Some mathe- 
matical details are relegated to Appendices. 



II. QUALITATIVE DISCUSSION 

The purpose of this Section is to describe qualitatively 
how the the non-equilibrium fluctuations affect the den- 
sity of states of the superconductor. We need first to 
classify the collective excitations, which are present in 
the system, and to understand how non-equilibrium con- 
ditions influence them. 

We start by recalling the basic physics of phase fluctu- 
ations in superconductors, then consider their coupling 
to currents in the metallic layer. The electric current in 
the normal metal is accompanied by shot-noise. As we 
show below, this noise gives rise to the classical phase 
fluctuations. Finally, we demonstrate that the enhanced 
fluctuations dramatically affect the BCS density of states. 
In particular, they lead to a non-zero density inside the 
BCS gap. This Section is concluded by suggesting an 
independent experiment to test our theory. 



A. Collective modes in N-S sandwich 

Consider a superconducting film at zero temperature. 
It is well known, that all of the excitations with the en- 
ergy smaller than the superconducting gap A are asso- 
ciated with the phase 9 of the order parameter. The 
time evolution of this phase is governed by hydrody- 
namic equations, which in the absence of external mag- 
netic fields can be written asH 



3 s = 0, 
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(2.1a) 



(2.1b) 



(2.1c) 



where n s is the perturbation of the carrier density in 
the superconductor, v s is the thermodynamic density 
of states per unit area in the s uper conductor, and j a 
is the supercurrent. Equation (2.1a) is th e con tinuity 
relation. We wrote the London equation ( 2.11 ) for a 
dirty superconductor and expressed the superfluid den- 
sity through the diffusion coefficient D s in the normal 



state of the superconductor. Equation (2.1c) is the con- 
ventional Josephson relation between the electrochemical 
potential and the phase of the order parameter 9. Finally, 
the electrostatic potential tp is connected to the density 
variation by the Coulomb law 



e<p= dr'V{r -r')n s {r'), V(r) 



(2.2) 



Performing the Fourier transform of Eqs. (2.1) and 
( |2.2D , we obtain the dispersion relation for the collective 
mode 



= T AD S Q 2 [1 + v s V{Q)] , V(Q) = 
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It corresponds to the usual 2D plasmon with dispersion 
u> ~ \fQ which has little effect on the behavior of the 
system because of its small density of states at low fre- 
quencies. 

When a layer of the normal metal is brought nearby, 
V(Q) gets screened, and the dispersion relation Eq. (2.2) 
becomes linear. To see this, one has to include the nor- 
mal cur rent s into the set of the hydrodynamic equat ion s 
( |2.lD - (2.2). The equation for the scalar potential (2.2) 
is modified to 



etp= / dr'V(r-r')[n s (r') + n m {r')] 



(2.4) 



dispersion. This is due to the fact that the currents in the 
normal metal flow along the interface in direction oppo- 
site to the currents in the superconducting layer, so the 
total charge accumulation does not occur. The physics of 
this mode is essentially similar to the well known Schmid 
- Schoixj mode in the vicinity of the critical temperature 
or to the second sound in supcrfmidal The only dif- 
ference is that the normal excitations are not thermally 
activated in the superconductor itself but rather exist in 
the normal metallic layer close to the superconductor. 
This, however, has little consequence on the charge dy- 
namics. 



Here we neglected for simplicity the thickness of the iso- 
lating layer between the normal metal and the super- 
conductor assuming that d <C 1/Q. A ll the results are 
insensitive to this assumption (see Sec. |III C ). The den- 
sity of carriers n m in the normal metal is governed by 
the continuity equation and the Ohm's law: 
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(2.5a) 



(2.5b) 



The charge in the normal metal is redistributed 
by the electric field according to Eqs. (2.5). We 
evaluate n m fr om these equations, substitute the re- 
sult into Eq. (2.4), and find, that the potential be- 
comes dynamically screened. At frequencies much 
smaller, than the plasmon frequency in the normal layer, 
uj p ~ v m Vo(Q)D m Q 2 , the screened potential takes the 
form 



V(Q,lj) = 



1 —iuj + D m Q 2 



(2.6) 



We now substitute Eq. (|2.6|) into the dispersion law of 



the collective mode Eq. (2.3) and find 



Jph 



u ph 



Q 



I, _ * 



ttAA 
h 



1/2 



A 

n' 



Vrr, 



1/2 



(2.7a) 

(2.7b) 
(2.7c) 



Equations fl2.7(fl - ( |2.7b| ) are valid provided that u) 1 



ph 



> 



u P h- 



This condition is satisfied already at small frequen- 
cies Tito ~ hiv' ph ~ A(G S /G m ) <C A. The lifetime of this 
mode is finite due to the interaction with the relaxation 
mode in the normal metal. 



Since 



uJ 



ph 



ph ' 



Eqs. (2.7b) and (2.7c) describe a 



well pronounced collective m ode. We will call this mode 
"phason". According to Eq. (2.7b), the phason has linear 



B. Phase fluctuations due to the current in the 
normal layer 

We have seen that a presence of the metal gives rise to 
the new collective mode, the phason. This mode corre- 
sponds to the oscillating electric currents flowing in the 
opposite directions in the metallic and the superconduct- 
ing layers. 

Now let us consider what happens, when a current is 
driven in the normal layer. The average currents in the 
metal are accompanied by the fluctuations known as the 
shot noise. Since the currents in the metal are coupled to 
those in the superconductor, it is natural to expect that 
in the superconductor the fluctuating currents appear as 
well, and consequently, the phasons are generated. In 
other words, a (ic-current in the metallic layer should en- 
hance phase fluctuations in the superconductor. 

To include these fluctuations in our description of the 
N-S sandwich we add Langevin sources 5j l to the current 
in the normal metal. Equation (2.5t) takes the form 



j m = -<J m Vtp - DmVri m + Sj[ 



(2.8) 



The fluctuations 5j l are described by their correlator 
(&jf$jl)u>,Q- Provided the frequency u> is much less than 
the applied voltage u <C eU/K and the energy-relaxation 
is negligible, this correlator can be written asS 



{f>j?&3i)u,Q = Sai3<J m eU cx e(j r 



(2.9) 



The superconducting phase 9 in the presence of the 
current fluctuation s 5j , ca n be det ermi ned fro m th e sys- 
tem of equations (2.1), (f2~4|), ( 2.5a ) and (|2^). In 
addition, we use the Einstein relation er m = e 2 v m D m . 
As a result, we can present the phase fluctuation 89 as 
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(2.10) 



where uj 2 h (Q) is the phason dispersion 



In Eq. ( |2.iq ) 

and all the subsequent formulas, we have not specified 
an inessential numerical prefactor, which will be found 
in the next section. Using the correlator Eq. (2.9) we 
obtain 
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Therefore, the correlator of the phase fluctuations 
{6 2 ) u q has a well pronounced phason pole and is pro- 
portional to the applied voltage U. 

In what follows, we will need the single point corr elato r 
of the phase fluctuations, i.e., the integral of Eq. (2.11) 
over the momentum Q. The main contribution to this 
integral comes from the pole, which corresponds to the 
resonant excitation of the phason by the current fluctua- 
tions in the n ormal layer. The logarithmic divergence at 
Q = in Eq. (2.11) is not important and, as we show in 
Sec. IIIB, disappears due to the gauge invariance. Since 



the linewidth of the phason decreases w ith t he increase 
of the diffusion coefficient D m , see Eq. (2/7c), the large 
factor D m in the denominator of Eq. (2.11) is canceled. 
As a result, the single point correlator of the phase fluc- 
tuations (SO 2 )^ does not depend on the parameters of the 
metallic layer 



(SO 2 ), 
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(2.12) 



Equation ( [2.12; ) is valid, provided G S A/G m <C twj < eU, 
where G Sjm denote dimensionless conductances of the su- 
perconducting (in the normal state) and normal layers 
respectively: G sm = 2itTio s 



Je 2 = 2irTiv^ m D^ m , these 



are the conductivities measured in units of e /2nh 
1/(25.8/07). 



C. Effect of the phase fluctuations on the tunneling 
DOS of the superconductor 

We have found that in the presence of the normal 
layer the phase fluctuations in the superconductor are 
large due to the resonant excitation of the phasons. Now 
we are interested in the effect of these fluctuations on a 
measurable quantity, e.g. on the tunneling conductance, 
Gt = 01/ dV of the junction, with I and V being the 
current and the voltage across the junction respectively. 

In the lowest order in the tunneling amplitude, the 
tunneling conductance Gt is determined by the density 
of states of the superconductor Eq. (1.2). The density 
p(e) depends on single particle excitation energies. In 
the absence of fluctuations the excitation energy in the 
superconductor is given by the usual BCS expression 
E = \/ + A 2 , where £ is the energy of the orbital state 
counted from the Fermi level. The energy E can not be 
smaller than A. This prevents electrons (or holes) from 
the metal with energies smaller than A from tunneling 
into the superconductor. 

However, in the presence of the phase fluctuations, it 
becomes possible for an electron with the energy smaller 
than A to tunnel and then to absorb a phason to com- 
pensate for the energy deficit. As a result, the density of 
states turns out to be finite even inside the gap E < A. 



To describe this effect of the phason assisted tunnel- 
ing more quantitatively, we first calculate the density of 
states in the superconductor in the presence of homoge- 
neous phase fluctuations. In this case, due to the orthog- 
onality of the orbital wave functions, all the transitions 
are confined to the same orbital (characterized by some 
orbital energy £). The problem simplifies since within the 
single orbital we have to consider only four states (one 
orbital can be occupied by no more than two electrons): 
tp - empty orbital, ^2 filled orbital, Vt and _ singlets. 
States tp and ip 2 are coupled to each other due to the 
exchange with the condensate, whereas singlets are not. 
The resulting Schrodinger equations are 



ihii> = Ae 2l9 (^ 2 , 



ihip 2 = 2£02 + Ae 
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^0, 



(2.13a) 



(2.13b) 



(2.13c) 



At frequencies smaller than ell/ft, the occupation num- 
ber of phasons is large ~ eU/Twj, and, therefore, 6 can 
be treated as a classical variable. 

Phase 9 changes with the characteristic fre quen cy of 
the order of eU/h <§C A/h. Hence, Eqs. ( 2.13 ) can 
be solved in the adiabatic approximation. The time- 
dependent ground state of this four-level system is given 
by 



*GS(*) 



Wo + vip 2 e 2l6{t) 



-j- J* dtiEaaiti)^ (2.14) 



where u and v are the usual coherence factors, u 
v = sin a and tan 2a = — A/£,and 



£ = t-ne(t). 



cos a, 



(2.15) 



So far we discussed an isolated superconductor. Let 
us now consider tunneling of an electron with the energy 
e from the normal metal into the superconductor. Since 
the initial state has the energy Egs + 6 and the final 
state is one of the states V'TU)' each having energy £, the 
transition amplitude of this process can be estimated as 



A(t) = 




Eg sit") - § 



(2.16) 



The prefactor in Eq. (2.16) denoted by . . . includes the 
tunneling matrix element and the coherence factors. This 
prefactor is not important for present discussion. The 
transition rate is given by 



-J-oc hm \(\A{t)\ 2 ) 



(2.17) 
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where (. . .)g stands for the averaging over the fluctuations 
of the phase 9. 



We can substitute the amplitude Eq. ( 2.16 ) into the 
tunneling rate Eq. (2.17) and use the explicit form of 
Eas(t) given by Eql~f 2.15). In the exponent we have 



The density of states in the superconductor p(e) is pro- 
portional to the tunneling rate l/r tun „. Evaluating the 
time integral in the exponent, we obtain 



p(e) = lim — ( 

t^oo ZTTt 



(2.18) 



The absolute value square of the time integral can be ex- 
pressed in terms of the double integral over t 1 and t" . It is 
convenient to change variables to the sum and difference 
of t' ± t" . The exponent depends only on the time differ- 
ence. Indeed, the averaging over fluctuations of 6 can be 
performed independently and the average functions de- 
pend on the time difference only. The integral over t' + 1" 
then yields t, which cancels with the denominator in the 
prefactor. In the remaining integral over the difference 
t' — t" the upper limit can now be taken to infinity. The 
last step is to sum over all orbitals, since the electron can 
tunnel to any of them. This amounts to integration over 
£ — y v s f d£)- This integral is independent from the 

orb 

fluctuations of 6, since they produce shift of £ which is 
irrelevant because of the integration over £, and yields 
the BCS density of states po{e) Fourier transformed to 
the time domain: 



-A3)« 
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Thus, we find 
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Therefore, in the presence of the phase fluctuations the 
density of states in the superconductor becomes modified 
by the factor, which describes the phason assisted tun- 
neling from the normal layer. To ev aluate the average 
we use the phase correlator Eq. ( 2.12 ). It is not correct, 
since here we considered only the homoge neou s fluctua- 
tions of the phase, while the correlator Eq. ( 2.1 2| ) includes 
also inhomogeneous fluctuations. However, the effect of 
the inhomogeneity can be estimated as HD S Q 2 / A, while 
the main contribution comes from the phason mode with 
D S Q 2 ~ hu; 2 /A, Therefore the correction is ~ ?i 2 u 2 / A 2 , 
and should be neglected because the corrections of the 
same order were already omitted within the adiabatic 
approximation. 



Averaging the phase factor with the help of Eq. (2.12) 
we find 



oo oo 



oo 



exp 



eU 
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(2.20) 



The integration is cut off at to — eU, because at to > eU 
the classical description fails, and Eq. ( 2.12| ) for the phase 
fluctuations is not valid. 

In the energy interval A — eU < e < A one can expand 
the exponent to the first order, which corresponds to a 
single phason process. As a result, 



eU 
2G S A 



(2.21) 



where numerical coefficient is the result of rig orous treat- 
ment of the following section. Equation ( 2.21 ) is the main 
qualitative result of the paper. It shows that in the pres- 
ence of the normal layer driven out of equilibrium by 
the applied voltage U, the superconductor no longer has 
a hard gap. Instead, there is a dip, and the tunneling 
density of states is suppressed at e < A (see Fig. p. 



P(E) 




A-eU 



FIG. 1. Sketch of the tunneling density of states p(e) the 
superconductor modified by the current fluctuations in the 
normal metal driven out of equilibrium by the applied volt- 
age U. 

At energies smaller than A — eU the density of states 
differs from zero due to mul ti-ph ason processes. It means 
that the exponent in Eq. (2.20) should be expanded to 
higher orders. However, to describe these multi-phason 
processes, we need better understanding of the nature of 
the frequency cut off. Also, we have so far neglected the 
quantum part of the phase fluctuations, to > eU. These 
issues are addressed in the technical part of the paper, 
Sec. Ill where we perform a rigorous calculation of the 



tunneling, based on Keldysh diagrammatic technique. 
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D. Is the effect of phason assisted tunneling 
independently observable? 

We have shown that the phason assisted tunneling 
manifests itself in the tail in the density of states within 
superconducting gap. When interpreting the experimen- 
tal data, this effect could be confused with the broaden- 
ing of the distribution function /m(e) in the normal metal 
[see Eq. (1.2)]. We believe that such a misinterpretation 
of the experimental data was made in R ef. B. 

We should warn the reader that Eq. (2.2 1|) as well as 



the subsequen t rigo rous calculation of the tunneling cur- 
rent [see Eq. (|3~6g| ) for the result] were obtained for the 
simplest model of the N-S junction, namely the 2D sand- 
wich. The spectrum of collective modes is sensitive to 
the geometry of the system, therefore, our results are not 
expected to describe the experiment of Ref. || in detail. 
However, we have presented a strong evidence of the ex- 
istence of the microscopic mechanism responsible for the 
change of the shape of the tunneling conductance, which 
is different from the trivial broadening of the distribution 
function in the normal metal. In this paper we do not 
intend to evaluate the effect of phason assisted tunneling 
in various possible experimental realizations of the N-S 
junction. Instead, we suggest an independent measure- 
ment, that can distinguish between the two mechanisms. 




FIG. 2. Scheme of the experiment for observation the effect 
of non-equilibrium phase fluctuations on the density of states. 
Switch "SI" corresponds to the experiment of Ref. ^[ Switch 
"S2" corresponds to the measurement of the density of states 
affected by the current fluctuations in non-equilibrium elec- 
trode "Ml" but not convoluted with its distribution function. 

The suggested experimental setup is shown on Fig |^. 
This scheme differs from the experiment of Ref. |l| by 
adding another normal electrode "M2". The electrode 
"Ml" is driven out of equilibrium by the applied voltage 
U. That affects the density of states in the superconduc- 
tor ("SC") due to phase fluctuations as well as broadens 
the distribution function in the metal "Ml" itself due to 
energy relaxation. We suggest to measure the tunneling 
conductance between the superconductor and the second 
electrode "M2" . Since this electrode is in equilibrium, its 
distribution function is broadened by temperature only. 
Such measurement will only show the modification of 



the density of states of the superconductor by the cur- 
rent fluctuations in "Ml" (shown on Fig. |l|). Comparing 
the two measurements, one should be able to determine 
which effect dominates the tunneling between the super- 
conductor and the non-equilibrium normal metal. 



III. THEORY OF THE INTERACTION EFFECTS 
IN N-S JUNCTIONS. 

In this Section we set fi = 1 in all intermediate expres- 
sions. 



A. Tunneling current 



In this section we use KeldyshE£JEil diagrammatic tech- 
nique to evaluate the tunneling current I between the 
superconductor and the normal metal. In order to take 
superconductivity into account properly, we need Green 
functions Q to be matrices in the space K ® N, which 
is the dkect product of Keldysh space K and Gor'kov- 
NambuEij space N. We choose the basis in the Keldysh 
space, for which electronic Green's functions aretU 



9 = 1 Q A 



(3.1) 



K 



Let [ ]_ and [ ] + denote a commutator and anticom- 
mutator correspondingly. The three components of the 
Green's function Q a p, where a and /3 can be either M 
(metal) or S (superconductor), are defined as 



05j(l, 2) = -i v (t! - ta)<[*a(l) ® *j,(2) 
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(3.2c) 



where we used a sho rt-hand notation n = (f n ,t n ) for 
n = 1,2. In_Eqs. (|3.2|) T](x) is Heaviside step function, W 



are Namb spmors 



* = 



and ip^(n)i ip-\(ri) are the fermionic operators in the 
Heisenberg representation. Since we are using the 
Keldysh technique, (. . .) denotes the average with arbi- 
trary distribution function. 

We start our calculation of the tunneling current I be- 
tween the superconductor and the normal metal with the 
general expression 



-T Q Tr[a 2 {g M s-QsM)]; 



(3.4) 
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Within our choice of the basis [Eq. 3.1 the current vertex 
in Eq. is Pauli matrix a 2 in Keldysh space 
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(3.5) 



TV in Eq. ( p.4[ ) means trace in the K ® N space and 
also assumes that r\ = r 2 = r, where r corresponds 
to the contact point. Below we often ommit the spa- 
tial coordinates of the Green functions and write explic- 
itly only temporal coordinates. It always assumes that 
r i = r : 2 = f. 

In the first order in the tunneling amplitude To the 
Green's functions can be calculated as 

g MS = -ir 2 v m v s T J dt' SM (t,t')e- ieVt ' S s(t',t), (3.6) 

where V is the applied voltage and gM and gs are the 
Green's functions in the normal metal and in the super- 
conductor, respectively, defined as 



la(tl,t 2 ) = Gaa(f,r;t 1 ,t 2 ) 



(3.7) 



In what follows the semiclassical Green's functions in the 
K®N space will be denoted by boldfaced symbols, while 
the 2x2 matrices in the Nambu space will be denoted 
by a hat. 

In the normal metal the semiclassical Green's function 
gM satisfies the Usadel equation 

- D m V c R (g M o V c r Sm) 

where the matrix r 3 is 

1 0\ „ fl 



T 3 = 



o D K -\o -i /A 



(3.9) 



the commutator designates 

[*,gJ|f]t = *(tl)gM(tl,t2)-gJtf(tl,*2)*(t2), (3-10) 

and 

gM °gM(ti,t 2 ) = J dt 3 g M (t 1 ,t 3 )g M (t 3 ,t 2 ). (3.11) 
The scalar potential is also a matrix 



ip+ if— 





i 




J ® 

' K 








(3.12) 



TV 



While the uniform scalar potential can be gauged out 
and doesjipt affect the behaviour of the system, it can 
be shownEj that the slow spatial variation of ip can be 
taken into acount by means of linear functional K[tp] so 
that the resulting Green's function can be written as 



iK[<p]{ti)r 3 



(gMO + gMl) e" 



-iK[ v ](t 2 )r 3 



(3.13) 



where gMo(^i — t 2 ) is the noninteracting Green's func- 
tion, while gMi(ti,t 2 ) contains all of the corrections not 
captured by the transformation Eq. ( 3.1 3| ) . 

The functional K[<p] can be specifically chosen in such a 
way, that in equilibrium the correction gMi(ii, t 2 ) is por- 
portional to the square of the gradient gMi ~ {^K[cp]) 2 . 
When the system is out of equilibrium, gMi(^i,^2) also 
acquires a term proportional to the deviation of the dis- 
tribution function /m(c) from equlibrium one fM,eq(e), 
i-e. SAfir^ [(p](f M (e) - f M ,eq{e))- Such functional is 
given byEj (see also Appendix A) 



K = 



K[<P]+ 
K[<p]- 

D m Q 2 ) 




= K 



V- 



jy(u>) 



-{iLU + D m Q 2 Y x 



(3.14a) 



(3.14b) 



where N(uj), which can be named "bosonic distribution 
function" is defined as 

N(u) = [ - [2/ M (e)/ M (e + w) - f M (e) 

J U! L 



-/m(£ + w) 



(3.15) 



In this case the corrections gMi(ti,t 2 ) in Eq. ( 3.13j ) 
come from the effects of energy relaxation and interac- 
tion localization corrections (the effects of weak localiza- 
tion were negleced from the very beginning, as discussed 
above). These corrections are inversely proportional to 
the metal conductance G m (in 2D; in ID the corrections 
are ~ l/G m (l p h) where l p h is the dephasing length; for 
details see Ref. 4). Therefore, these corrections can also 
be neglected. 

The retarded and advanced Green functions g^o m 
the time domain are delta functions 



4(<) = ~9Mo(t) = *(*)■ 



(3.16) 



The Keldysh function in the energy domain is related to 
the distribution function by 



ffM0 (e) = 2(1 - 2/ M (e)). 



(3.17) 



Fluctuating electric fields give rise to the fluctuations 
of the phase 9 of the order parameter A in the supercon- 
ductor . The fluctuations of the absolute value of the or- 
der parameter A do not couple to the fluctuations of the 
electric field, and their spectrum has a gap ~ A. There- 
fore they can be ignored. The superconducting Green's 
function gs can be evaluated by solving the Usadel equa- 
tion in the superconductor 
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where 



D S V C R (SS o V c flgs ) 

Ae ie * 



r 3^5 + ^ r 3„ i[(A + $)!gs]t = 0i (3 _ 18) 



A = 



-Ae~ ieK 



N 



and 



'A' 



(3.19) 



(3.20) 



A 



In Appendix B we show that at small frequencies 
(ui ~ V <C A) the dominant effect is captured by the 
gauge transformation 



gs(*l,*2) 



o i01t(ti)T ! 



gso(*i - h)e 



-i8 K (t 2 )T 3 



(3.21) 



where gso is the usual BCS Green's function. 

We now substitute the Green's functions in the nor- 



mal and sup erconducting layers Eqs. ( 3.13 ) and ( 3.21 ) 
into Eq. (3.6) and substitute the re sult into the expres- 
sion for the tunneling current Eq. (3.4). The next step 
is to average over the phase and Coulomb fluctuations. 
The current is then given by the integral 



/ = 



Go 



dt'e- ieVt 'Tr 



{<J 2 [gMo(t-t')g s (t' -t) 



-Es(t - t% M0 (t' -t)]}, 



(3.22) 



where G = 2iru m v s e 2 TQ. We used the cyclic property of 
the trace and an obvious fact that the matrices <72 and 
t 3 to obtain the averaged function 

gs(ti - h) = (e^^gsoCfi ~ h)e-^ T3 ), (3.23) 
where the fluctuation field equals to 

cf) = 9 - K[tp}. (3.24) 



So far our analysis was quite general and independent 
of the geometry of the sample. Even though we were 
consid ering tunneling through point contacts, Eqs. ( 3.22| ) 
and (3.23) can be strightforwardly generalized for the 
wide contacts with large number of channels. 

We now calculate the trace in Keldysh space and use 
t he e xplicit for m of the metallic Green's functions Eq. 
(3.16) and Eq. (3.17) to evaluate the time integral. The 
differential tunneling conductance Gt = dl/dV can be 
written as 



Gn 



Gq 

~2~ 



de df. 



M 



2tt de 



(e-e^)ReTr[gf(e)r 3 ]. (3.25) 



The last step is to perform the average in Eq. (3.23). 
In the leading approximation in 1/G m the fluctations are 



Gaussian, and we obtain (see Appendix C for details) for 
the function hf(e) = Tr[gf (e)r 3 ] 

OO 

h§(e) = J dte™ [/i+-(t)e ai ( D w 





■A^WeHVO-VI")] , (3.26) 



where 



h+-{t > 0) = -ho + (-t) = —K^iAt) (3.27) 

and K\ is the modified Bessel function. The fluctuation 
propagators are defined as 



M>+) = ^(«i-*2), 

i+^-) = «Z&(*i-*2), 
!»-</»-) =0, 



and 



V 



+-(-+) 



(3.28a) 
(3.28b) 
(3.28c) 
(3.28d) 

(3.29) 



The conductance can finall y be expressed in terms of 
the function h§(e) Eq. ( 3.26 ) as 



G f de 8/m i inn tr, \ 



(3.30) 



So far our results are quite general. They describe 
the effect of the low frequency (lj < A) collective mode 
on the tunneling conductance. In the follo wing sections 
we evaluate the tunneling conductance Eq. (3.3C) for our 
particular setup. To do this we need the fluctuation prop- 
agators T>, which is evaluated in the following subsection. 



B. Propagators for low energy excitations. 

Throughout this paper we are using the small param- 
eter 1/G m . When the fluctuations propagators T>^ is 
calculated, this parameter allows us to restrict ourselves 
by ladder diagrams, i.e., to use RPA. These diagrams 
can be summed up by means of solving the Dyson-type 
equation 



V = V n + VnlW. 



(3.31) 



In this equation the propagators T> and the polarization 
operator II are 8x8 matrices defined in the K®MS®tpd 
space. The 2x2 space <p9 describes the two fluctuation 
fields we have in the probl em. Since these fluctuations 
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take place in normal metal as well as in the superconduc- 
tor, another 2x2 space MS is needed. 

In the metal sector only the ipip polarization operator 
is present, which yields the usual diffusion pole 



(n 



MM\R 
<P<P ' 



D m Q 2 ~ i~ 



(3.32) 



where D m is the diffusion constant in the normal metal. 

In the superconductor there exist the phase as well as 
the Coulomb fluctuations. Corresponding polarization 
operators can be obtained from the general relation (om- 
miting the SS superscript) 



-U^ik = tt^Tt [rf g s ] + 5 iv 5 ktp , 
v s bipl 



(3.33) 



where the indices v can be + or — . The polarization 

ik 



operator UP? is the matrix in Keldysh space: 



n 



ik 



(3.34) 



Indices i, k (which can be ip or 9) label the fluc tuation 
parameter space. The second term in Eq. ( p. 33 ) is the 
"anomalous" contribution of the large e region which is 
not captured by the Usadel equation (for details see, for 
example, Ref. 5). The vertex Ti is a vector in the ip8 
space (as indicated by index i) and a matrix in K ® N 
space given by 



r£ = iA(r") 



K 





,-2i0 



„2i6»r. 



N 



r£ = i{r") K ® In, 



(3.35a) 
(3.35b) 



where (t+) k = 1 K and (r ) K = {<t 2 )k- 

To proceed with this progr am and evaluate the polar- 
ization operators Eq. (3.32) we need an explicit form 
of th e Gr een's function gs- The approximate solution 
Eq. ( 3.23 ) which was obtained by the gauge transforma- 
tion is not sufficient here, since the polarization operator 
is a gauge invariant quantity. For this reason substitu- 
tion of Eq. ( p. 23] ) into Eq. ( [3.33j ) immediately gives 
zero. Therefore we have to take into account the gradi- 
ent terms, which were neglected in Eq. ( 3.23 ). To do 
that we expand the Usadel equation Eq. (3.18) to the 
first order in fluctuation fields 



£ ) s Q 2 gso(ei)gsi(ei,e 2 ) 



-«[Ho,gsi] e = i[6U,gso] 



(3.36) 



where H = Ik <E> H , SU = # + 2iA® + SA. The last 
term is irrelevant, since the fluctuations of the absolute 
value of the order parameter are gapped. Nambu space 
matrix Hq is given by 



e A 

-A -e 



(3.37) 



The gapless phase fluctuations have the following matrix 
structure (see also Appendix B) 













' K 




i) 



(3.38) 



N 



The solution of the first order equation Eq. ( |3.36 ) is 
the following Keldysh matrix (omitting the subscript S 
herefrom) 



(3.39) 



The anomalous function g x is only present before aver- 
aging over the fluctuations. 



where 



B_+ = D S Q 2 - iS-(ei) - iS+{e 2 ), 



S±(u) = \/(e± iO) 2 - A 2 . 
The rest of the functions are given by 



(3.40) 

(3.41a) 
(3.41b) 



9?- 



B++ 



+gf tanh^(5 + ( £l )-^( £l )) 



(3.42) 



D- 



g$6H + g£ - 5H+ + g£5H.g« 



(2 



+g? tanh^(5 + ( e2 )-5_( £2 )) 



(3.43) 



and 



g£SH + g£ +g£5H + g£ 



+g£5H_g£ + g«6H_g* - SH. 



:,A" 



+^tanh^(S + (e 1 )- 1 9_(e 1 )) 



(3.44) 



A' 



i n the obvious notation (compare with Eqs. (3.40) and 
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Using the first order sol ution we can calculate the po- 
larization operators Eq. ( |3.3S ) for the phase and the 
Coulomb fluctuations in the superconductor. For the re- 
tarded component of the polarization operator we obtain 



,2 r 



SS\R 



(IT") 



vs 



1 —iu) 

ILU LU 2 — nAD s Q 2 



(3.45) 



The solution of Dyson equation Eq. ( 3.31 ) is 

©=[©o 1 -n]" 1 , (3.46) 
where II is the polrization operator of the form 



n 



n ss 
o n MM 



(3.47) 



MS 



The unperturbed propagator T>q has non zero elements 
only in the Coulomb interaction channel 



V = -V(Q) 



1 

-Qd 



e -Qd 
1 



1 



K 



SM 



1 




(3.48) 



where the Coulomb interaction in two dimensions 
V{Q) — 2-Ke 2 /Q and d is the thickness of the sample. 
This form of the unperturbed propagatoris model depen- 
dent. It corresponds to the particular setup of the S-N 
sandwich, which we are discussing in this paper. 

To obtain the retarded propagators we have to invert 
the retarded block of t he 8 x 8 matrix. Since we need only 
<p<p propagator in Eq. (3.26), there are three elements of 
the remaining 4x4 matrix which are relevant: 



Vff = — 
M 



i + U vv ! IT 



n^n, 



V 



MS 



1 

M 



1 - Qd 



(3.49a) 
(3.49b) 



W = — 

s M 



V 

(1 - Qd) 2 



■~ + n 



M 



V 2 



(3.49c) 



where 



n 



M 



i + vn M 



V = 4ire 2 d. 



(3.50a) 
(3.50b) 



Equations (3.49) are written for the case Qd <C 1. 

Sub stituting t he po larization operators Eqs. (3.45) and 
(|3 32[) into Eqs. fl3.49|), we obtain 



Tiff — 
U M — 



M 



- ttAD s Q 2 1 



1 



V 



MS 



1 iujv s 



qjee 1 \_ QD m Q 2 



IUJV S 



M V -iu + D m Q 2 ' 



(3.51a) 
(3.51b) 

(3.51c) 



where £ = v n 
form 



v s {\ + v m V) and Eq. (3.50a) takes the 

D m Q 2 



V *-iw + D m Q 2 



P(Q 2 ) = Q 2 



P{Q% 



(D m t:A(D s 



(3.52a) 
(3.52b) 



The pole P(Q 2 ) in the fluctuation propagators means 
that there is a well-defined collective mode with linear 
dispers ion r elation. This mode is the phas on, di scussed 
in Sec IIA. The imaginary term in Eq. (3.52t) deter- 
mines the finite lifetime of the phason. However, this 
term is small, since D m /D s ^> 1, and the phason life- 
time is large by this parameter. 

We now combine these propagators into a single re- 
tarded propagator of the field <p using the definition of 



(Kl<p}K[c P ])-2(K[cp]8) + (68). 



(3.53) 



Substituting the explicit form Eq. ( 3.14 ) of the lin- 
ear functional of the Coulomb fluctuations in the normal 



metal and the propagators Eq. ( 3.51 ) we find 



v 



v s nAD s C P(Q 2 ) 



V V 
ttAD, 



-iui + D m Q 2 



D m Q 2 



(3.54) 



Note, that although each of the propagators Eq. (3.51) 
contains the 1/Q 2 factor, the combined propagator Eq. 
([D34j) does not. This is a manifestation of the gauge in- 
variance: the zero mode here corresponds to a uniform 
perturbation with the constant scalar potential, which 
can be gauged away and can not affect the physics. 

Finally, to obtain the time-dependent pr opaga tor, 
which is neede d to calculate the average Eq. ( |3.26 ) we 
integrate Eq. ( 3.54 ) over the momentum. The last two 
terms in brackets do not contribute by their analytic 
properties, while the first term gives for the retarded 
propagator 



2AG, 



The corresponding Keldysh propagator is 



ISgTLUJ 

AG S 



N(co), 



(3.55) 



(3.56) 
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where the dimensionless conductance of the superconduc- 
tor G s = 2ttu s D s . Th e dis tribution function N{u>) can 
be obtained from Eq. (3.15) using the particular form of 
the metallic distribution function /m(c) corresponding to 
the case when the system is driven out of equilibrium by 
the applied voltage U 



f M (e) 



eU 



eU 



7) e- 



This results in the bosonic distribution function 

if M > eU 



ujN(lu) 



UeU+\w\) , if \u\ < eU 



(3.57) 



(3.58) 



of intera ctio ns. Then it is the BSC Green's function given 
by Eq. (CI), which at frequencies near the gap diverges 
as ~ 1/ y/\A — if. The equilibrium correction Eq. ( 3.63| ) 
at long times behav es lik e l/(At) (since the upper inte- 
gration limit in Eq. ( |3.63 ) is cut off by the gap). Expand- 
ing the exponential will only give rise to extra factors of 
| A — e| in the numerator, thus being small. Therefore we 
can neglect the equilibrium correction. 



The non-equilibrium contribution in Eq. (3.26 ) is t hus 
dominant and the average Green's function Eq. ( 3.25 ) at 
e > is 



Reh§ 



eU eU 



e - A ' 7rAG. 



(3.64) 



C. Results for tunneling conductance. 



Using the propagators Eqs. ( 3.55 ) and ( 3.56 ) we now 
calculate the average Eq. ( 3.26 ). First, let us separate the 
equilibrium (V = 0) and non-equilibrium contributions 



(3.59) 



The non-equlibrium part comes entirely from the 
Keldysh propagator Eq. ( 3.56] ): 



eU 



K 



eU=0 



Taking the Fourier transform we get 



V ne {t) 



ieU 
2ttAG s 



Si(eUt) 



sin eUt 
eUt 



- 1 



(3.60) 



(3.61) 



where S-Uf) = lnz+7 — Ci(z), Ci(z) is the integral cosine 
functionE_3 and 7 = 0.577 ... is the Euler constant. 

The remaining equilibrium part is determined by the 
Fourier transform of 



T>eq(t) 



1 r 
2 



*W) - 2?^(0) 

which is given by the integral 

00 

dio 



1 



eU=0 



V$Jt), (3.62) 



V eq {t) 



AG S J 2tt 




1 



1 



2ttG, G.Ar 



(3.63) 



Here we have used the fact that the advanced propagator 
vanishes at t > 0. 

It is easy to see that when substitued in to Eq . ( [3.26 ), 
only the non-equlibrium propagator Eq. (3.61) gives a 
noticable deviation from the non- i ntera cting (Golden- 
rule- type) behaviour. Consider Eq. (3.26|) in the absence 



where F is a function of two dimensionless variables 
y = eUj (e — A) and z — eU '/irAG s . It is given by 



F(y,z) = Rc 




e ixsgny e -z[S 1 (\y\x) + ^ 



(3.65) 



The tunneling conductance Eq. ( |3.30|) d epends on the 
real part of the Green's function Eq. (3.64). The external 
voltage in the experimental setting is always less than the 
superconducting gap, so that z«l. To probe the vicin- 
ity of the gap we need to look at e ~ A corresponding to 
the limit y 1. In that limit Eq. ( |3.65|) gives 



F(y;z) = \y\ 



f , if y < 
1 



if y > ' 



\y\ » 1. (3.66) 



That describes the behaviour of the density of states 
in the superconductor in the vicinity of the gap edge, in 
particular the appearance of non-zero density of states 
inside the gap (y < result) 

For completeness we include the solution in the oppo- 
site limit i/< 1 which describes the tail of the density of 
states as we go deeper into the gap. When e > A we can 
expand the exponential to obtain 



F(y;z) 



1 + T& zy\ 



y>0; y<l. 



(3.67) 



However, deep in the gap the perturbation theory gives 
zero (due to energy conservation : no single particle 
process is possible at e < A — eU). Therefore, we 
have to evaluate the full integral. This yields the non- 
perturbative result, which describes the multi-phason 
processes 



\y\z 1 



(3.68) 



Here e = 2.718. . .. 

Finally, we ex press the result for the tunneling conduc- 
tance Eq. ( |3.30| ) in terms of the function F 
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A 



=eV±eU/2 



eU eU 



lei — A' 7rAG. 



(3.69) 



where F is given by Eq. ( 3.65|) and it s asymptotic behav- 
ior is described by Eqs. (p!66Q-(|3.68D. 
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APPENDIX A: 



IV. CONCLUSIONS 

In this paper we have studied the effect of electron- 
electron interactions on the tunneling conductance of N-S 
junction with the metallic layer driven out of equilibrium. 
Contrary to the common belief, the tunneling current can 
not be presented as a convolution of bare BCS density of 
states and the derivative of the non-equlibrium distribu- 
tion function in metal. We have demonstrated, that the 
interaction between the tunneling electrons and those in 
the normal metal drastically affect this current modifying 
the superconducting density of states. This modification 
is manifested in particular in appearance of finite d ensity 
of states at energies even smaller than A, see Sec. Ill C . 
This effect can complicate the clear determination of the 
distribution function by tunneling experiments similar to 
Ref. § 

The reason for the strong effect of the non-equlibrium 
state in the metal on the superconducting density of 
states is apperance of low-energy collective mode on the 
N-S junction, analogous to the Schmidt-Schon mo de or 
the second sound, which we called "phason", Sec [I A. 



The tunneling can be accompanied by emission or ab- 
sorbtion of phasons. The shot noise in the normal layer 
generates these phasons, thus, affecting the tunneling 
density of states. 



The particular form of the phason spectrum Eq. (2.7) 
and conse quent ly the result for the tunneling conduc- 
tance Eq. ( [3.69] ) depend on the geometry of the junction. 
In this paper we considered the simplest model of the 
junction, namely the 2D sandwich. Calculations of the 
effect of the phason assisted tunneling on various exper- 
imental realizations of the N-S junction are outside of 
the scope of this paper. Therefore our, results can not 
be expected to describe the experiment!] in detail. How- 
ever, we have demonstrated an existence of a microscopic 
mechanism changing the tunneling conductance, which is 
different from the trivial broadening of the distribution 
function in the normal metal. To test which mechanism 
dominates in reality we have suggested the independent 
experiment, see Sec. II D 
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Here we demonstra te th e choice Eq. ( 3.14 ) of the func- 
tional K[ip] in Eq. ( |3.13|). Wc simply substitute the 
Gree n's function Eq. ( |3 . 13 ) into the Usadel equation 
Eq. ( |3.8| ) and require that in equilibrium the correction 
SMi(ti,t2) is porportional to the square of the gradi- 
ent gMi ~ (V-RT[<p]) 2 and out of equilibrium to the vari- 
ation of the distribution function gA/i ~ AK[tp](fM (e) — 

/M,eg(e)). 

To proceed with this program we start with the trans- 
formation 



Sm 



(Al) 



After the substitution into the Usadel equation Eq. (3.S) 
we find the equation for the function h.(ti, t%) 



<9h dh 



dt 



dh dh 
*,h] t =0, 



(A2) 



where the covariant derivative now includes the commu- 
tator of the gradient of K[<p] and the function h 



V c R h = V R h + i[V R K[<p],h]. 



(A3) 



Here we assumed that the functional K[tp] (which is a 
matrix in the Keldysh space) commutes with the scalar 
potential matrix <fr which we will verify later using the 
explicit form of K [ip] . 

We now treat Eq. ( |A2| ) in perturbation theory h = 
gA/o + gMi, wher e gjy/ o is the free Green's function de- 
termined by Eq. (3.16). The equation becomes 



D m / dt 3 g M o(h - t 3 ) V 2 K[cp}(t 3 )T 3 S M0(t 3 - t 2 ) 



dK[ip] 
dh 



-*(ti) g M o(*i-*a) 



T 



gun</*./>)('^-*(t 2 ; 



(A4) 



The functional T in the right hand side of Eq. (A4) 
(which explicit form is not needed for what follows) con- 
tains higher powers of K[<p], We now choose such a form 
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of K[ip) that all linear terms (the left hand side of Eq. 



(A4)) cancel with the external potential To do that 
we treat the left hand side of Eq. (A4) as the equation 
for K[ip] (with zero in the right hand side). 

To simplify the solution we notice that the non- 
interacting Green's function has the matrix structure 



gMO = 9mo 



(A5) 



where guo is the 2x2 matrix in Keldysh space and r 3 acts 
in Nambu space. Therefore each term in the equation is 
proportional to r 3 which therefore can be omitted. Then, 
using the normalization gjvfo ° gi/o = 1 an d performing 
the Fourier transform in space variable (only K[ip] de- 
pends on i?), we obtain the equation for the functional 
Klip] 



dt 3 gMo(h - t 3 )K[(p](t 3 )g M0 (t 3 - t 2 ) 

-K[ip]{tx)8{t x -t2] 
( dK[ip] 



- <p(ti) gMo(h - t 2 ) 



„ , , fdKlip] „, . 
-9M0[h - h) \ —gf 2 y(*2. 

where ip is the matrix 
and the functional K[ip] has the same structure 



= 0, (A6) 



(A7) 



K[if] 



K + K_ 
K- AT 4 



K 



The Green's function gMo is given by 
9Mo(h-h)-\ Q - S ( tl -t 2 ) 



(A8) 



(A9) 



A" 



where the Keldysh co mpone nt is related to the distribu- 
tion function (see Eq. ( 3.17 ) ). 



We substitute Eqs. (A7) - (A9) into the matrix equa- 



tion Eq. ( |A6| ) and investigate each component separately. 
The diagonal elements ("11" and "22") of Eq. @) differ 
only in sign and give the equation for A'_ 



-g^-iP-(t) 



0. 



(A10) 



which gives the solution in frequency domain (see Eq. 

(B ) 



AT_ 



<P- 



iuj + D m Q 2 



(All) 



This solution also s atisfi es the lower non-diagonal 
("21") element of Eq. (A£). The upper non-diagonal 



(Keldysh) element gives the equation for K + which we 
write in Fourier space 

(D m Q 2 - iu)K+{u)[f{e) - f(e - u)] 

+2D m Q 2 K_(Lu)ip_ [/(e)/(e - w) - 1] 

= V + [/(e)-/(e-w)]. (A12) 
The last trick is to introduce the bosonic "distribution 



funct ion" N(lo) Eq. ( 3.15 ). In the second term of Eq. 
( |A12| ) we rewrite 

f(e)f(e - w) - 1 = -N(u)[f(e) - f(e - w)] 



+ [N(u;)[f(e) - f(e - w)] + f(e)f(e - w) - lj . (A13) 

In equilibrium the expression in square brackets is 
equal to zero, which is just a reflection of the de tailed 
balance principle. Then, using the solution Eq. (All) 
for AT_, we obtain the solution 



K+ = 



9^ 



-iw + D m Q- 



- 2^_Re 



N(u) 



ioj + D m Q 2 



(A14) 



which was previously given in the matrix form in Eq. 



Out of equilibrium the solution Eq. ( A14 ) is not exact, 
since the detailed balance princ iple is no longer valid and 
the square bracket in Eq. ( A13 ) gives rise to corrections, 
which are proportional to the difference between equilib- 
rium and non-equilibrium distribution functions. These 
corrections should then be included in the next order 
Green's function gMi- 

The calculation of the correction gMi introduces the 
Altshuler-Aroonov corrections to the conductity as well 
as the energy relaxation. We will not dwell on this issue 
in the present paper. 



APPENDIX B: 



Here we show that the gauge transformation Eq. (3.21) 
captures the dominant effect of phase and Coulomb fluc- 
tuations in the superconductor. Since we integrate the 
fluctuation propagators over momentum before using us- 
ing them to calculate the tunneling conductance, we can 
here set the momentum to zero from the very b eginn ing 
taking the integrated Keldysh propagator Eq. (3.56) as 
known. We only need the Keldysh propagator since we 
are interested in non-equilibrium effects (see the para- 
graph following Eq. ( p. 63] ). In the straightforward per- 
turbation theory described in Sec. IIIB| the terms de- 
scribing the phase fluctuations contained the large factor 
of A. These contributions can be taken into acco unt b y 
means of the gauge transformation (similar to Eq. ( |3.2l| )) 
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gs(*i,*2) = e 



g(ti,te)e 



-iG K (t 2 )T 



(Bl) 



where now the function g(ti,i2) is not the BCS Green's 
function but is to be determined from the Usadel equa- 
tion, which after the transformation becomes 



D s V c (g o V c g) + [H , g] e + [6H, g] e = 0, 



(B2) 



where 



m, S ], 



de 3 
2ir 



<5H(ei - e 3 )g(e3,£2) 



-g(ei,e3)<JH(e 3 - e 2 ) 



(B3) 



the covariant derivative after the gauge transformation 
includes a commutator with the gradient of 9 



V c g = Vg + i[V6» A -r 3 ,g] 



(B4) 



and <5H = $ + 0. 

In addition to the equation Eq. (B2) the function 
g(ii,t2) has to satisfy the constraint 



(B5) 



which fixes the normalization. 

Now instead of A the phase fluctuation term has a 
factor of u>. In what follows we show that all the correc- 
tions calculated by perturbation theory will be small at 
least as a power of uj/A and thus the gauge transforma- 
tion indeed captures the dominant contribution. Since 
we already have the factor of w n the numerator, the 
only possibility to obtain a strong correction is to ex- 
cite a soft mode. Therefore to simplify the discussion we 
first separate the gapped modes which can not give rise 
to strong corrections (since we are interested in energies 
much smaller than A) and then treat the remaining soft 
modes in perturbation theory. 

The separation of the modes can already be seen in the 
absence of fluctuations. The free (BCS) solution is given 
by Eq. @ and can be written in the following way 



e<73 + Acr 2 
v /(e±«0) 2 -P 



(e±z0) - A 



(0"3 + 0"2 ) 



-(03 - (72 ) 



(B6) 



where only in this section we denote by 02 and 03 the 
following matrices in the Nambu space 



c 2 





-1 



0-3 



A' 





-1 



(B7) 



A 



Here the large and small parts of the solution (in the 
limit e — A <C A we are interested in) turned out to have 
different matrix structure. 



To see this structure in the equation we first explicitly 
separate the fast oscillating part of the Green's functions 
writing the solution as 



(B8) 



Here we are focusing on positive energies. 

Now the functions gi varies slowly in time (since the 
external potential SH is a slow function) and we get the 
equation 

- £> s V c (gi o V c gl ) + [H, gl ] £ + 0H, gl ] £ = 0, (B9) 

where the Hamiltonian H has an additional term 

H = H +iAl K ®cr3- (BIO) 



We now see that the large term in the BCS Green's 
function Eq. (B6) commutes with the Hamiltonian Eq. 
( B10| ) (neglecting the e term for the moment), while the 
small term anticommutes. In the same manner we will 



look for a solution of Eq. (B9) which contains matrices 
commuting and anticommuting with H 



gi = Gi ® ljv + G 2 <S> (03 + o- 2 ) 



+Hi (g) a\ + H 2 (g) (cr 3 - cr 2 ) 



(Bll) 



where hat denotes matrices in Keldysh space, while the 
a j act in Nambu space. The anticommuting terms Hj 
are small and the commuting terms Gj are large. In the 
absence of fluctuations G\ = Hi = and Gi and H2 



form the BCS Green's function Eq. (B6) 



When we substitute the solution Eq. (Bll) into the 
transformed Usadel equation Eq. (B9), in the first order 
we only use the non-commuting terms Hj in the commu- 
tator with the Hamiltonian, neglicting the smaller contri- 
bution from the other terms in the equation. As a result 
we obtain 

(e 1 -e 2 )[G 1 <g> cr 3 + Gi <g> 1 N ] + (ei + e 2 )G 2 ® a x 
+[«5H, Gi] e ® 1 N + [5H, G 2 ]e O (ff 3 + ^2) 



+2AH ± ® (cr 3 + o- 2 ) - 4A^ 2 ® <J\ 



= 0, (B12) 



where • • • denote gradient terms which we do not include 
here due to their complexity but will restore in the next 
equation. 

Collecting terms with the m atrix structure correspond- 
ing to four matrices in Eq. ( Bll ) we obtain the final 
equations for Hj and Gj 
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(ei - e 2 )Gi - D S VA 4 



-D S [V9 K ,(A 1 +A 2 )} = 0, 



(B13a) 



A 



2 A y (e ± id) ' 



(B16c) 



(ei - e 2 )G 2 + [SU,Gx} e - D S VA 1 - iD s [V9 K , (A 4 + A 3 



2D s [Ve K ,G 2 oVe K G 2 }=0, 



(B13b) 



[8U,G 2 ] e + 2AH 1 - % -D a [V6 K ,Ax] - ^D S {V6 K ,A 2 } 



D S V(A 3 + 2%G 2 o V0 K G 2 ) = 0, 



(ei + e 2 )G 2 - AAH 2 - D S VA 2 + iD s {V9 K , A 4 } 



(B13c) 



iD a [V9 K , (A 3 + 2iG 2 o Ve K G 2 )} = 0. (B13d) 



Here the curly brackets denote anticommutators. The 
matrices Aj are the gradient terms porportional to G\ 

A x = Gi o VGi + iG x o [V9 K , G 2 ] 



iG 2 o[Ve K ,Gx] 



(B14a) 



A 2 = -iGx o [V9 K , G 2 ] - iG 2 o {V6 K ,Gx}, (B14b) 



A 3 = Gx o VG 2 + G 2 o VGx + A 4 , 



A A = -Gxo[V6 K ,Gx]. 



(B14c) 



(B14d) 



Since all the terms in Eq. (B13a) are proportional to 
Gx it is still not generated in the first order, therefore all 
the Aj terms in the first order are equal to zero. The 
only correction from the gradients comes from the term 
quadratic in G 2 , which changes the equation for G 2 in 
the presence of phase fluctuations 

(ei - e 2 )G 2 - 2[X79 K , G 2 o X79 K G 2 ] = 0. (B15) 

However this is the equilibrium correction: for the only 
term coupling to the non-equilibrium distribution func- 
tion, the classical field 9 + , the commutator is equal to 
zero. In any case, at zero frequency the correction is 
equal zero and thus it does not produce any additional 
singularity. Thus the first order corrections are small by 
a factor of 1/A 



Gi =0, 



G 



R,A 



(e±ioy 



(B16a) 



(B16b) 



Hx = -^[(SH, G 2 } ( + 2iD s G 2 o V 2 9 K G 2 , (B16d) 

wher e e is counted from A after the transformation 
Eq. ©. 

Thus we see that 5H does not couple soft modes to 
each other and so it can introduce only perturbative cor- 
rections of the order (<5H/A) 2 . Therefore in the leading 
order in 1/A the dominant effect of the phase fluctu- 
atio ns is indeed captured by the gauge transformation 
Eq. (Hfl) . 



APPENDIX C: 



Here we p erfor m the averaging over the fluctuating 
field in Eq. ( 3.23 ). The BCS Green's function gso has 
the following structure in frequency domain 



H 



9so ( £ ) 



Sff (e) = tanh -^-( 5 f - g£ Q ) 



± iO) 2 — A 2 

e 

2T' 



(Cla) 
(Clb) 



where the matrix in Nambu space is 
Ho 



e A 
-A -e 



(C2) 



sc 



In time domain the matrix structure is the same and we 
can represent gso as a sum of two parts 



9so = h r 3 + 1 t 2 



(C3) 



To calculate the tunneling current we only need the nor- 
mal component of gso, therefore we can focus on averag- 
ing ho- Then the matrix t 3 can be carried through, so 
that 



h = (e^h Q (t 1 ~t 2 )e^ t ^). 



(C4) 



Since only the retarded function enters the expression 
for the tunneling current Eq. ( 3.3C| ) , we now multiply the 
Keldysh space matrices to get h^. The Green's function 
has it's usual (rotated) form 



ho 



h R h 
n n 



K 



h J K 



while the fluctuating field is 



(C5) 



(C6) 



K 
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After matrix multiplication we get for the retarded func- 
tion 

hfiih-h) = (eJ^+^^+^lh^ cos ^(h) cos <p^(t 2 ) 



h cos</>_(ii)isin</>_(£ 2 ) 



+Hq sin <j>_ (ti) sin0_(t 2 ) 



(C7) 



By definition a retarded Green's function is only non- 
zero when its time argument is positive. In that region 
an advanced Gre en's function is zero. Therefore we can 
simplify Eq. ( |C7| ) as 

h$(h > ia) - h e i (*+W-++(*>'» cos<M*i) 

\h§ + hg)e-*+- (t2) + (h$ - h$)e 1 *- ( * 2) 1 . (C8) 
We are now ready to perform the average: 
/ e *(<M*i)-<M*2))+i(</>- (ti)-<t>- (t 2 ))\ _ 

= e (<P+<P+)-(<t>\)(0) + (4>-<P+) + (<P+<P-) ^ ^ C9a ) 
/ e *(^+(tl)-*+(*2))-»(^-(ti)+*-(*a))\ _ 

= e <0+0+)-«>(O)-(0-0 + ) + (0 + 0-) j ( C9b ) 



where all the correlators (4>4>) depend on the time differ- 
ence t\ — t 2 , except where indicated. 

We now introduce propagators T> of the fluctuating 
fields 



b + <f> + )=W^{t x -t 2 ), 
b + <j>J) = W% (t x - t 2 ), 



(ClOa) 



(ClOb) 



(ClOc) 



At ti > t 2 the advanced propagator is zero. There- 
fore only sums and differences of th e re tarded and the 
Keldysh functions are present in Eq. flC8| ). These can be 
expressed in terms of functions of the original Keldysh 
basis as tiff ± h$ = 2h^ ' + -* and similarly with the 
fluctuation propagators to obtain the final expression in 
time domain 

h§(h > t 2 )= h+-{h - t 2 )e 2i ( v t;^-^-v:*W) 

-ho + (ti - kJe^V^-^-Vw) . (Cll) 

When Fourier transforming to th e freq uency domain, 
we notice that the first term in Eq. ( Gllj ) contributes at 
positive frequencies, while the second term contributes 
only at negative frequencies. Therefore 



ft*(e > 0) = / dte let h+-(t)e 
o 



ut h+-(t)e 2i ( v ^^- v ^^) . 



(C12) 



For reference purposes we list here the explicit form of 



the BCS function ht in the time domain 



h+-{t > 0) = — Ki(iAt), 
where K\ is the modified Bessel function. 



(C13) 
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